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04 . ABSTRACT 

We address the problem of angular momentum transport in stellar radiative interiors 
py^ with a novel semi-analytic spectral technique, using an eigenfunction series expansion, 

^jy that can be used to derive benchmark solutions in hydromagnetic regimes with very 

high Reynolds number (10 7 — 10 8 ) . The error arising from the truncation of the series is 
evaluated analytically. The main simplifying assumptions are the neglect of meridional 
circulation and of non-axisymmetric magnetic fields. The advantages of our approach 
are shown by applying it to a spin-down model for a 1 M Q main-sequence star. The 
^ . evolution of the coupling between core and envelope is investigated for different values 

\ of the viscosity and different geometries and values of the poloidal field. We confirm 

that a viscosity enhancement by ~ 10 4 with respect to the molecular value is required 
to attain a rigid rotation in the core of the Sun within its present age. We suggest 
i that a quadrupolar poloidal field may explain the short coupling time-scale needed 

OO ' to model the observed rotational evolution of fast rotators on the ZAMS, while a 

i dipolar geometry is indicated in the case of slow rotators. Our novel semi-analytic 

*~0 ' spectral method provides a conceptually simple and rigorous treatment of a classic 

MHD problem and allows us to explore the influence of various parameters on the 
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| rotational history of radiative interiors. 

0^ ! Key words: MHD - methods: analytical - methods: numerical - stars: rotation 

stars: magnetic fields - stars: late-type 



1 INTRODUCTION 

The rotational evolution of a solar type star during the main sequence lifetime is deeply affected by the level of internal 
coupling that is established between the radiative core and the convective envelope at various ages, which competes with the 
angular momentum loss produced by a magnetised stellar wind. In spite of the observational constraints gathered so far, the 
identification of the physical process(es) ensuring a uniform rotation of the core, i.e., a substantial rotational coupling by the 
age of the Sun, is still lacking. 

Observational information to constrain the rotational evolution of low mass stars currently comes from two major sources: 
helioseismology and rotation period surveys in open clusters. The former provides evidence of a latitudinal dependence of 
solar rotation in the convection zone (CZ), the existence of a thin shear layer (tachoc line) at the top of the radiative core, 
and an almost uniform rotation at greater depths (at least down to 0.2 — 0.3 Rq, see Thompson et al. 2003). On the other 



hand, rotation period measurements of open cluster members effectively constrain the rotational evolution because of our 
fairly reliable estimate of their age. A satisfactory theoretical explanation of these data is still missing owing to the limited 
knowledge of the different processes involved in the evolution of stellar angular momentum. 

During the pre- main sequence (PMS) contraction phase, stars are s pun up, but c ontemporarily suffer a remarkable angular 



momentum draining by the interaction with their circumstellar discs (|Konig]Hl99lD . Very young (viz., a few M yr old) open 
clust ers already show significant spread and mass-dependent structures in their rotation period distr ibutions jLamm et aT 



20041). On the main sequence (MS), a ngular mome ntum loss by a magnetised wind in solar- type stars ijWeber &^D*£wijjl^)67r 



Kawalei 1988h was early recognised by Schatzman ( 1962 ) as the main cause of their rotational evolution, contrasting with the 



more massive stars' behaviour (see KraftJI 19671 ; Skumanich 19721 ) 



Surface period measurements are sensitive to the internal rotational profile. Phenomenological models of MS angular 
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momentum evolution ( Allain 19981 ; Bouvier 20081 ) have been successful in reproducing the observations if a certain amount 
of differential rotation between the radiative core and the convective envelope is allowed. In those models the exchange of 
angular momentum between the core and the envelope is governed by a single parameter, i.e., a coupling time-scale r c which 
fixes the rate at which angular momentum is transferred between the two regions to establish rigid rotation. The observed 
distribution of the rotation periods in open clusters can be reproduced by assuming that r c is of the order of ~ 10 Myr for 
the stars that begun their evolution on the ZAMS as fast rotators (i.e., with an initial rotation period of a few days), while 
it is remarkably longer (r ~ 100 Myr) for slow rotators (i.e., with initial periods of several days). An understanding of the 

processes eventually ensuring rotational coupling between core and envelope and within the core itself is still lacking; 

The issue is also closely related to a number of currently open questions in stellar evolution theory. For instance, 



Bouvier 



( 2008h tentatively explained the observed correlation between Li over-depletion and the presence of hot Jupiters around 
MS stars on the basis of these rotational evolution models because an initial slow rotation is a characteristic of stars with 
exoplanets, owing to their prolonged interaction with a proto-planetary disc during their PMS phase. Such objects are 
characterized by a long coupling time (r c ~ 100 Myr) which leads to the development of a sizeable differential rotation at the 
core-envelope interface. A sizeable shear produces an enhancement of the turbulent mixing which in turn makes Li destruction 
more efficient. To address the internal coupling processes from a theoretical point of view, ISpruitl (|l999l , l2002h discussed the 
stability of toroidal magnetic fields in radiative stellar interiors. He suggested that, due to the vertically stabilizing effect of 
the subadiabatic stratification of the core, motions associated with hydro-magnetic instabilities are constrained to the locally 
horizontal direction and their mean effect can be described through an enhancement of the momentum and magnetic field 
diff usivities up to turbulent valu e s seve ral orders of magn itude greater than their molecular counterparts. 



Charbonneau fc MacGregor ( 19931 ). hereafter CM93I . solved the equations for angular momentum transport and toroidal 



magnetic field evolution in a solar-type star throughout the MS lifespan, assuming axial symmetry and a pre-existent (fossil) 
poloidal magnetic field. They used a finite elements technique, finding that a sizeable toroidal field develops from the winding- 
up of the poloidal field by the differential rotation. When it becomes strong enough, the Maxwell stresses begin to react to 
any further amplification of the field and the system enters into a new regime characterized by torsional Alfven oscillations. 
This phase displays a remarkable energy dissipation as soon as the oscillations along neighbour magnetic field lines get in 
opposition of phase (phase-mixing), ow ing to their s li ghtly d ifferent periods, as a consequence of the density stratification and 
poloidal field gradient inside the core l|Spruijll999h . lCM93l find that the poloidal field geometry is an important feature to 
determine the time-scale for angular momentum redistribution inside the core, because phase mixing (much more efficiently 
than diffusion) enhances angular momentum redistribution and toroidal field reconnection on surfaces of constant poloidal 
field, thus producin g a quasi-stationar y regime with an angular velocity almost constant on them, according to Ferraro 
isorotation theorem (|Mestel et al. 19881 ). Another key finding is the impo ssibility to achieve a state of u ni form ro tation of the 
core within the solar age with molecular viscosity alone. Later work by Riidiger fc Kitchatinov 1 199rJ ) I RK96h . based on a 
finite difference numerical approach, confirmed these results and remarkably the "viscosity de ficit" p roble m. 

We solve the same coupled, non- homogeneous partial differential equations (PDEs) of CM93 and RK 9~rj for angular 
momentum tran sport and toroidal field evolution, with the technique of eigenfunction expansion ijMorse fc Feshbachl Il953l : 



Habermanl 12004 ) . We find an exact analytic solution of the problem that can be expressed as a series expansion. Formal 



substitution of the series into the equations yields an infinite system of linear, first order ordinary differential equations (ODEs) 
that can be truncated according to the required degree of accuracy. To this purpose, we complement our implementation with 
an analytic formula for the truncation error. In principle, our spectral method has no limitation in terms of Reynolds number 
and can be used as a benchmark to test the accuracy of other numerical methods to solve the same problem. This new rigorous 
treatment of a classic MHD prob lem is t he n ovelty reported in the present work. We shall compare the results of our approach 
with those of previous works bv lCM93l and lRK9d briefly discussing its advantages for the study of the angular momentum 
transport in a radiative stellar core. Our treatment is particularly relevant for an accurate study of the phase mixing process 
for which previous numerical techniques were not capable of a rigorous treatment of the contribution of subgrid lengthscales. 



2 MODEL 

We model the MS evolution of specific angular momentum in the core of a 1 Mq star solving the angular momentum transport 
equation and the toroidal component of the induction equation. 



2.1 Basic assumptions 

We consider an inertial reference frame with the origin O at the barycentre of the star and the polar axis z along the rotation 
axis. Spherical polar coordinates are adopted, r being the distance from O, $ the colatitude measured from the North pole, 
and ip the azimuthal angle. 

The present calculation is based on a number of simplifying assumptions. We restrict our domain to the radiative core, 
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assuming that the angular momentum transport time-scale in the CZ is much shorter than in the radiative interior, as 
turbulent viscosity there is at least 10 — 12 orders of magnitude greater than molecular. The system is assumed to be strictly 

axisym metric. 

The lSpitzerl (|l962h expressions for molecular viscosity and magnetic diffusivity are: 

= 1.2 ■ 10" 16 T 5/2 



"mol 



_ i 9 _ 1 

"p cm s 



77 mol = 10 13 T- 3 / 2 cm 2 s" 1 . (1) 

As pointed out by RK961 molecular viscosity is far too low to reconcile the differential rotation regime that is likely 
established in young, solar-type stars with the unifor m rotati on of the core of the present Sun as deduced by helioseismology. 
To overcome this "viscosity deficit" problem, following RK96I . we introduce an artificial viscosity enhancement parameter Red 
such that: 

V c ff — Ref£V mo h (2) 

This suffices to the purpose of pre senting our method and comparing the results with previous works. ISpruitl |l999h and 
iDenissenkov fc Pinsonneaultl l|2007h present a possible identification and physical explanation for such an effect. 

Application of Eqs. fl} requires the knowledge of the stellar structu re (viz. the depth dependenc e of de nsity and tem- 
perature). Here we use the so-called model S of the Sun introduced by Ichristensen-Dalseaard et al. 1 19961 ) for the whole 
computation. The evolution of the stellar structure on the MS has a minor impact on our computations and can be safely 
ignored. 

We neglect the Eddington-Sweet circulation because of the extremely long time-scale that is of the order of 10 12 yr. Note 
that a differential rotation which is not uniform along cylindrical surfaces around the z axis drives a meridional circulation 
because of the non-potential character of the associated centrifugal force. We shall neglect such a circulation because the 
subadiabatic stratification of the core effectively opposes motions in the radial direction, strongly reducing its velocity. 

Our assumptions rule out the possibility of poloidal field regeneration by dynamo action a nd 3D magnet ic instabilities. 
According to the analysis of magnetohydrodynamic (MHD) instabilities in radiative regions by Spruit (1999), however, the 
fastest-growing instability should be the m = 1 mode of the Tayler instability, which may be included as an additional Maxwell 
stress term in our equations. 

Our assumption of an axisymmetric poloidal magnetic field is justified by the fact that in a differentially rotating core any 
initially non-axisymmetric field component is smoothed o ut by winding up and diffusion on a time-scale r w much shorter than 
the diffusion time-scale of the axisymmetric component. ISpruitl (| 19991 . Sect. 3.1) estimates r w of the order of 10 2 yr for the 
Sun on the ZAMS. Note that if the initial field is so strong to oppose the w inding up by the initial differential rotation, it can 
possibly find a non-axisymmetric equilibrium (|Braithwaite fc SpruitJl2004r ). However, in this case the core would be rotating 
rigidly from the outset being completely coupled by the strong field and the angular momentum would be transported on the 
Alfven time scale which is of the order of 10 2 — 10 4 yr, i.e., much shorter than the wind braking time scale. We do not treat 
this case because such a strong coupling inside a young solar-like star is not in agreement with the phenomenological models 
of rotational evolution requiring coupling time scales of 10 7 — 10 8 yr to account for the observations, as discussed in Sect. [T] 



2.2 Governing equations 

With the hypotheses discussed above, the total fluid velocity u and magnetic field B can be written as: 
u = r sin $Q(r, t)e v , 



B = B„ + Bt = 



r sin i? 



13f(r,i),{) <9*(r,tf,t) 



dr 



-et) 



+ B v (r, t)e v . 



The axisymmetric magnetic field is expressed in terms of two scalar functions, the flux function ^ and the toroidal component 
B<f,. The poloidal field lines lie on surfaces of constant <& (magnetic surfaces) as it is immediately apparent from B p ■ VP = 0. 
The evolution of Q, is governed by the angular momentum conservation law (Rudiger 1989h : 



pr 2 sin 2 i9^p = V- (pr 2 sin 2 feVfi) + ^-B p ■ V(r sin&B v ), 

at 47r 

and that of B v and $ by the induction equation, which breaks into two scalar equations ( Radler 1980h : 



dB v 
dt 

a* 



v 2 - 



B v + 



V77 ■ V(r sini9B ¥ 



r 2 sin 2 i? 
a 2 * sintf d ( 1 <9* 



r sin •& 



+ rsintfBp ■ Vfi, 



(3) 



(4) 



(5) 



r 2 d~& ysini? cW 

The evolution of as given by Eq.j5j is decoupled from Eqs. [3] and [4] but it requires knowledge of initial conditions 
which are generally not available (cf. IRK96I) . 

The presence of a magnetic field within the core is expected as a relic of dynamo action during PMS. This problem has 
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been addressed by iKitchatinov et ajj (|2001f ) for stars like the Sun. Their numerical calculations go from an age of a few Myr 
to about 30 Myr and account for the dynamical retreat of the CZ during the PMS evolution. They find that non-axisymmetric 
modes are the most readily excited in young PMS models, but they are gradually replaced by axisymmetric ones as the star 
approaches the ZAMS. Moreover, any early non-axisymmetric field is rapidly diffused away, as already noted at the end of 
Sect. 12.11 We shall therefore consider a stationary, axisymmetric configuration for the poloidal seed field, assigned through a 
flux function ^, and specialize our considerations to the case in which $ can be factorized as ^(r, 8) = tp(r)x(9). Specifically, 
we investigate the case of a dipole confined within the spherical shell n ^ r ^ r% (see Sect. 12. 3[) with: 

* = BoRt (r - n) 2 (r - r 2 ) 2 sin 2 0, (6) 

where Bo is the scale of the field intensity and R* is the radius of the star. The next multipole, studied for comparison purpose, 
is the quadrupole, defined by: 

* = B Rl (r - r±f (r - r 2 f sin ■& cos ■&. (7) 



2.3 Boundary and initial conditions 



Thanks to the assumption of axisymmetry, our computational domain encompasses a meridional section of the radiative core: 
(r, 0) G [ri,r 2 ] x [0, 7r]; r 2 represents the lower boundary of the CZ (for the Sun, r 2 = 0.7 Rq), while n is a nonzero lower 
boundary, introduced to avoid singularity at the origin. The actual value of r\ does not significantly affect the results, provided 
that it is chosen in such a way that the cylinder where r sin 6 ^ n contains a negligible amount of angular momentum. We 
are interested in the global response of the core to an external wind braking torque, which we take into account through a 
suitable specification of the boundary conditions at r%. 

Specifically, on the radial boundaries we assign the angular momentum fluxes. At the inner boundary ri, we assume it 
to be negligible. At the outer boundary r%, the angular momentum flux is equal to that lost via the magnetised wind: 



Or 







dr 



Wit). 



(8) 



For simplicity, we neglect any dependence of the wind torque and dQ/dr on the latitude, which is in any case a second-order 
effect, given the high turbulent viscosity of the CZ which makes the amplitude of the latitudinal shear on top of the boundary 
significantly smaller than that in the core. For the specification of the function W(t), which includes a model of the wind 
braking process, see Appendix IA1 

For the toroidal magnetic field we use "insulating" boundary conditions: 



B- 







B- 



0. 



(9) 



They prevent the development of toroidal magnetic fields of unrealistically large intensities, as shown by, e.g., Garaud fc Guervillvl 
[ 2009T I m numerical simulations of the solar tachocline. It is interesting to note that tachocline models predict the existence of 
a circulation inside that layer that confines the interior magnetic field preventing its outward diffusion (e.g., ISpieeel fe Zahnl 
19921 : iGoueh fc Mclntvre!ll99"i : iGoughl l l2007f l. Nevertheless, even if the poloidal field diffused outward, the strong radial shear 
present in the tachocline would wind it up producing a strong azimuthal field. In the mildly subadiabatic environment of the 
tachocline, it would become unstable and emerge on a time scal e of 100 — 1000 days by, e.g., doubly diffusive instabilities 
i Schmitt fc Rosnerlll983l : ICaligari et Zlll995l : Isilvers et al1l2009al lbh. Therefore, such a toroidal flux will be rapidly removed 
from the upper boundary of our computational domain on a timescale so short as to justify the assumption that B v = at 
the outer boundary. 

Initial conditions for the problem at hand are the outcome of the PMS evolution, namely, the contraction of stellar 
radius from several to about one solar radii and the development of a convectively stable core. Although the angular velocity 
profile emerging from those processes is not known, we may assume that they establish an internal differential rotation with 
a gradient mainly in the radial direction. Therefore, we take as initial conditions: 



n in (r) 
fio 

Bin Jn 



l + 



0, 



11 ( r - r stcp 

— I atan I — - 

2 7r \ Ar s tc P 



(10) 
(11) 



where AQ is the amplitude of a smooth step in the initial rotation profile Qi n (r), centered at r s tep and with a width Ar s t ep 
(see Sect .[4j , Qo ~ 2.7- 10 -6 s _1 is the present solar angular velocity, used as a reference scale, and B Vt i n is the initial toroidal 
magnetic field. 

These assumptions for the initial conditions should be regarded as working hypotheses, merely chosen for the sake of 
simplicity. In any case, they have little effect on the solution, as it is the driving wind torque at the outer boundary r 2 that 
dominates the subsequent angular momentum evolution after a brief transient phase. 
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3 SOLUTION METHOD 

The method used to solve the coupled PDEs Q and Q is based on a novel, semi-analytic approach, relying on an eigenfunction 



expansion technique (see, for example. lHabermanll2004f ). Both Eqs. ((3j and Q can be put in the form: 



w(r,d)^u(r,d,t)-£[u(r,&,t)]=Q(r,#,t), (12) 

i.e. they consist of a linear homogeneous part (the l.h.s.), containing a differential operator C acting on the main dependent 
variable u, and of a r.h.s. source term (here called Q). C expresses the action of the spatial derivatives appearing in Eqs. 
([3| and (Q. Note that the source term in Eq. ((3]) depends on the assigned function W and the other dependent variable 
B v . Similarly, the source term in Eq. contains Q, thus producing a coupling between the two equations. Focusing on the 
homogeneous PDEs associated to Eqs. (|3} and (UJ), the technique of separation of variables is applicable to each of them. 
First, we determine the eigenvalues k„ and eigenfunctions <fi n (r , i?) of the linear differential operators C satisfying homogeneous 
boundary conditions, i.e., we solve the Sturm-Liouville problems determined by: 

C[<t> n ] + K n W 4> n = 0, (13) 

and ^f" = for Eq. <T3j , or 4> n = for Eq. Q on the radial boundaries ri, T2 and suitable regularity conditions on the z 
or 



axis. The function w that appears as a factor of the time derivative in the l.h.s. of Eq. (|12|l enters the eigenvalue equation 
(|13l) as a weighting function (in our specific case, w = pr 2 sin 2 i? for Eq. ((3]) and w = 1 for Eq. (0). If the differential operator 
C is self-adjoint (as it is the case for Eqs. © and Q, see Appendix [Bl for details), the solution u(r, $,£) can be represented 
as a series expansion in terms of the eigenfunctions <f) n , i.e.: 

oo 

u(r,&,t) = J2 a "(t)Mr,#), 

n=0 

where the coefficients a n (i) are functions of the time. In this way, the eigenfunctions provide the spatial dependence for the 
solution and warrant that it satisfies the homogeneous boundary conditions. 

The non- homogeneous source term Q, as well as the effects of non- homogeneous boundary conditions (as, in our case, 
the second of Eqs. ([HI), are taken into account through similar series expansions in terms of the same eigenfunctions because 
they form a complete set, as it is known from the theory of the Sturm-Liouville problem. In such a way, a system of linear 
ODEs for the evolution of the expansion coefficients a n (t) is derived by a direct substitution of the series expansions into the 
original complete equations or, more formally, by applying the Green identity, as it is explained in App.[B] This yields: 

^ +K„a„(t) = b n (t) + [...], 

where the terms in the r.h.s. indicated by the dots are those coming from the expansion of the non- homogeneous contributions, 
i.e., Q and the boundary conditions. 



3.1 Scaling of the equations 

Fixing a proper scaling for our variables is not a simple matter because of the remarkably different time-scales involved. Stellar 
rotation naturally sets the first time-scale: 

t n = 27rfi~\ 

where f2* is the mean surface angular velocity. Note that, for the whole MS, fl* is between 1 — 10 f^0, thus giving tn ~ 
10 -2 — 10 _1 yr. This is also the characteristic time of toroidal field winding up and amplification by differential rotation. The 
Alfven velocity ma introduces the Alfven crossing time-scale: 

B 

tA = -R./UA, 

where po is the value of the density at the top of the radiative zone (i.e., r = T2), chosen as a reference (numerically, tA ~ 10 3 
yr). In addition, global time-scales for momentum and magnetic field diffusion arise in our problem: 

t v = Rl/V0, 

tr, = Rl/rjo, 

with values of ~ 10 12 and ~ 10 10 yr, respectively, when the reference values of vo and r\o are taken at the top of the radiative 
zone. 

It should be stressed that tA, t v and t v are introduced here only for scaling purposes and are not representative of the 
time-scales of specific processes occurring in our system. 

The best compromise between such different time-scales is to choose tA as the unit of time and construct nondimensional 
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quantities as follows: t* — t/t\, r* — r/R*, B* — Bp/Bo, ^* — fiiA, B* — B v /B (cf. ICM93h . Dropping the asterisk 
superscripts on dimensionless quantities, the scaled equations become: 

1 



an 

pr 2 sin 2 tf— - -^-V ■ (pr 2 sin 2 tfz/Vfi) = B p ■ V(r sintfB^), 



at 



1 



v 2 - 



B v + 



r^dirBf 
r dr 



= rsini9Bp-Vn, 



(14) 
(15) 



where we exploited the fact that 77 is a function of r only, rj = drj/dr, and we introduced the following two "Reynolds 
numbers" : 
„ _ R*ua _ t v 

•^v — — 1 ) 

R*u A _ t„ 

^rj — ■ 
770 *A 



3.2 Separation of variables 



To perform the separation of variables in the homogeneous PDEs associated to Eqs. (|14|l and (|15p . we formally substitute 
factorised trial solutions: 

n(r,0,t) = w(t)Z(r,0), and 
B»,(r,0,t) = /9(t)S(r,t»). 
The standard separation procedure leads to: 

dtj A 
It 



— — 1- = 0, and 



V • (pr sin fivVZ) + Apr sin t?X 



0. 



(16) 



for Eq. (|14[l . introducing the separation constant A; and: 

1 " 



V ■ (t,V3) 



r 2 sin 2 i9 



0. 



(17) 



for Eq. (|15p . with p being the separation constant. Note that, since non-homogeneous source terms have been excluded for 
the moment, the solution asymptotically tends to f2 = const, and B v = 0, and the separation constants A and p, appearing 
as eigenvalues in Eqs. (|16[) and (|17[1 . have the physical meaning of reciprocal of the respective decay timescales. 

The spatial, two-dimensional eigenvalue problems can be solved by means of a further separation of variables, i.e., by 
assuming the following form for the Z and E eigenfunctions: 

B(r,0) = £(r)P*(0), 

where the angular factors P^ 1 ' 1 ^ and are Jacobi polynomials and associated Legendre functions, respectively, which satisfy 
the following equations: 



d_ 

dx 

d_ 

dx 



(1 



(1 



,dP, 



(1,1) 



dx 

dP 1 
dx 



+ 



+ n(n + 3)(1- x 2 )P, 



n(n + 1) - 



(1,1) 



P x 

1 n 



= 0, 
= 0, 



(18) 
(19) 



where x = cos 8, and n is a non-negative integer. Such angular eigenfunctions are the only solutions of Eqs. (|18[) and (|19|l 
that are regular at the boundaries x = ±1, i.e. on th e polar axis (for Jacobi polynomials and associated Legendre functions 



Smirnovll 19641 : lAbramowitz fc Stegunlll965l ) . 



see, e.g., 

After performing the separation of variables, the eigenvalues of the angular eigenfunctions, i.e., n{n + 3) and n(n + 1), 
appear in the respective radial equations with the role of parameters. In other words, to each value of n corresponds a radial 
eigenvalue problem for the functions Cnfc(r) and £„fc(r), as defined by the equations: 



r 4 dr 



pvr 



idC, n 



dr 



+ 



A„fcP - 



n(n + 3) 



pv 



1 d 


2 d£nk 






r 2 dr 


Vr dr 


+ 


pnk 



r 



n(n + 1) 



C,nk — 



0. 
0. 



(20) 
(21) 
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with the boundary conditions: 

Cnfctri) = Cnfcfa) =0, 

e»k(ri) = e«*(ra) =0, 1 j 

where the primes denote derivation with respect to r. The radial eigenfunctions have therefore two indices: n identifies 
the corresponding angular eigenf unction, while k = 0, 1, . . . labels the different radial eigenvalues and eigenfunctions in the 
respective sets (\ nk , ( nk ) and (fi nk , f nfc ). 

The radial equations (I20|l and (|21|l contain the functions p(r), i^(r) and r?(r) tabulated from the stellar model and must 
be solved numerically. The theory of Stiirm-Liouville problem ensures that the eigenvalues are real numbers, with a lower 
bound (i.e. a minimum eigenvalue that is A„o or fi n o) but no upper bound. The fe-th eigenfunction, belonging to the fc-th 
eigenvalue, is unique to within a multiplicative constant and has exactly k zeros in the open interval Jri,^ . 



3.3 Eigenfunction expansion 

For fixed n, the radial eigenfunction sets {£ n fc}, {Cnfe} are both orthogonal and complete in the interval (ri,r2), in the sense 
that any piecewise continuous function of r can be expanded as a generalised Fourier series of those functions in the same 
interval. The same is true of in and P^ with respect to their index n and for the angular variable x £ [—1,1]. 
It is therefore possible to write the general solutions of Eqs. (114p and (I15[) in the form: 

n(r,0,t) =J2 u >^(t)Uk(r)Pi 1 ' 1) (&), and 

n . k 

(23) 

B v (r,6,t) =X)/U(t)&*(r)Pi(0), 

n,fc 

respectively. The generalised Fourier coefficients are functions of time and can be expressed through the scalar products 
defined by the following integrals: 

-1 pr 2 

pr 2 sin 2 , &£ n kPi 1 ' 1 ^0,r 2 dr dcosi9 



u> nk (t) = &™L = JL^LI1 and (24) 



(Z,Z) 



2 

sm 



i 

+ 1 rr 2 



tilP^fdcosti / pr 2 {C, nk ] 2 r 2 dr 



j enfcW-PnW^(r,J?,t)r 2 dr dcostf 

P nk (t) = = J ^ Jr \. +1 ~ 2 , (25) 

J [P^f d cos d J [Ufr'dr 

where proper normalisation weights are introduced in the denominators, respectively. We recall that the radial eigenfunctions 
are orthogonal with respect to the weight function w(r) — pr 4 for Eq. (|20[) and w(r) = r 2 for Eq. (|21[1 . respectively, 
appearing in the corresponding eigenvalue equations. For simplicity of notation, in what follows we shall always assume that 
the eigenfunctions are normalised, i.e. (Z nk ,Z nk ) = 1 and {E. nk ,E. nk ) — 1. 



3.4 Complete time equations 

The equations determining uj nk {t) and f5 nk (t) as a function of t, including all the source terms in Eqs. (|14|l and (|15|l and the 
non- homogeneous boundary conditions, are derived in Appendix [B] making use of the Green formula for the bi-dimensional 
self-adjoint operators in Eqs. (|16[) and (|17[) . They read: 

" k + -^-LO nk {t) — Snkmh Pmh{t) +-^~^nk(t), 

at v v 

mh 

(26) 

^ mh 

This is a linear system of ODEs, containing 2 x N x K equations, where N and K are the number of angular and radial 
eigenfunctions retained in the series in Eq. (|23[1 . respectively. The system can be compactly rewritten as: 

x(()=Ax(i)+iw(!), (27) 

introducing the vector of the unknown functions, x(t) = (a>oo(t), • • • , ujmk [t), Poo{t), . . . , Pnk (t)) and the matrix A. Its prop- 
erties are discussed in Appendix [B] where we show that it is a normal matrix, i.e., AA T = A T A. For a normal matrix, a 



8 F. Spada et al. 



generalised spectral theorem holds, which states that it is always possible to diagonalise it by means of a unitary transforma- 
tion. The main difference with the well-known Hermitian case is that the diagonalised matrix is not necessarily real. In other 
words, the matrix A can be transformed into a diagonal matrix A through a linear transformation of the kind: A = U _1 AU, 
where U is a unitary matrix. Through such a transformation we reduce the solution of the set of Eqs. (|27|l to the integration 
of a decoupled set of ordinary differential equations of the kind: 

y i (t) = A ii y i (t) + g i (t), (28) 

with y = U _1 x, g(t) = U _1 w(f), and An the elements of the diagonal matrix A. Such equations are immediately integrated 
as: 

Viit) = yi{0) exp(-A«t) + f gi (t') exp[-A w (t - t')]dt'. (29) 

Jo 

Once the system (|28[) has been solved, we can go back to the original variables with the inverse transformation and 
therefore give the solution for the angular velocity and the toroidal field. 

As already noted in Sect. [3] the non- homogeneous source terms containing the poloidal magnetic field produce a coupling 
between u;'s and j3's, while the wind braking boundary condition has the effect of a driving force. Such a coupling allows 
the development of an oscillatory phase, despite the parabolic character of each of Eqs. (|I4|I and (|I5[) taken alone. From a 
mathematical point of view, the coupling terms give rise to the antisymmetric part of A that is responsible for the appearance 
of the imaginary parts of the eigenfrequencies leading to such oscillations. On the other hand, the symmetric part of A, 
which is diagonal, gives rise to the real part of the eigenfrequencies, which are negative, thus producing a damping of the 
oscillations themselves. Therefore, the method introduced makes the evolution of torsional Alfven waves more transparent 
from a mathematical point of view. Such oscillations are characterized by a periodic exchange of power between rotation 
and magnetic field, in the presence of a damping due to the diffusion of angular momentum and toroidal magnetic field (cf. 
Sect. E3| . 



3.5 Kinetic and magnetic energies 

The total rotational kinetic energy of the core is: 



Ikin - 2 



i [ pr 2 sin 2 titt 2 (r,$,t)d' i r = i f f pr 2 sin 2 d S~] ^nkC, n kPn' 1] V] UmhCmhPjn' 1 ' ] r 2 drd cos d 

ZJ ri J-i n>h m ,h 

= ^^2^2^" m (J P r2 (nk(,nhr 2 drJ io nk uj mh = i ^2 ^2S nm S kh uj nk uj mh = i ^2uj 2 k , (30) 

n.k m,h \ "1 ' n k m .h n,k 

where t> is the volume of the core and we have exploited the orthogonality properties of the angular and radial eigenfunctions, 
respectively (the sum is intended over n, k, both going from to oo). Similarly, the total angular momentum of the core at 
any time can be written as: 



Jcorc(i) = / pr 2 sin 2 m{r,$,t)d s r = J2\ f^^tiP^P^dcosti [ * pr A C, n kC,oodr 

„ t U-i Jn 



u) nk (t) _ woo(£) 



r p(i.i) r P (i,i)' 



where the orthogonality of the eigenfunctions was exploited and also the fact that £oo = const, and Pq = const. We see 



that the term involving ljoo in Eq. (|30[) is proportional to the total angular momentum, thus it is always nonzero if rotation 
is present. It is interesting to note that the state of minimum kinetic energy for a given total angular momentum corresponds 
to rigid rotation, with Tki n = f^oo- 

In the same way as the kinetic energy, magnetic energy associated to the toroidal field component can be calculated as: 

Ub = g^^[B v (r,^,t)] 2 d 3 r = -^/ftt- (31) 



n . k 



Note that Eqs. (|30ll and (I31p are clos ely related to the generalised Parseval relations for our Fourier series expansions 
:001; 



i Bovce fc DiPrimal koOI: IWeinbergerlll995h . that read 

/ pr sm ! m 2 {r,'d,t)drdcos$ = ^^Lo 2 nk {i), 

1 ''''l n=0fe = 

/ r 2 Bl{r,d,t)drdco S $ = ^^/? 2 fe (f). 

- 1 "'''1 n=0fe = 



The total energy contained in the system can be evaluated from Tkin and Ub, taking into account the proper scaling 
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coming from nondimensionalisation, as: 

where Eo = 5/0^0 an d -^o = poR* are the dimensional scales of energy and moment of inertia, respectively. 

In an isolated and ideal (i.e., with negligible dissipation) MHD system, the total energy and angular momentum are 
conserved. In a system with very high hydrodynamical and magnetic Reynolds numbers, energy conservation holds for a 
very long time, because dissipation is very small. We shall make use of this property to apply energy conservation to our 
system, obtaining a relationship for the truncation error in Sect. 14.11 Note that there is a minimum energy state in an 
isolated dissipative system which is set by the conservation of the angular momentum. It corresponds to rigid rotation with 
f2 = WqqCoo-Po and B v = 0. This is the final state attained by our system, when dissipation of the kinetic energy of 
differential rotation and magnetic energy of the azimuthal field is over. 



3.6 Numerical issues 



Our eigenfunction expansion technique requires the numerical calculation of eigenvalues and eigenfunctions of t wo Sturm- 
Liouville problems. This has been implemented using the sho oting method described b y, e.g., Press et al. 1 19921 ). For large 
values of n or k, it is possible to resort to the asymptotic theory ( Morse fe Feshbachlll953T l. The precision of these computations 
has been checked through the orthogonality of the eigenfunctions belonging to different eigenvalues, which is verified up to 10 
in relative units in our numerical integrations. One of the greatest advantages of our approach is that, once the determination 
of the set of eigenfunctions has been performed, the evolutionary calculations reduce to a matrix diagonalisation and inversion, 
and the integration of a set of de-coupled first order, linear ODEs whose general sol ution is given by Eq. (I29[) . Matrix inversion 
and diagonalisation are performed by means of the standard LAPACK routine^] ijAnderson et al. 19991 ). 

The possibility to estimate an upper bound for the error (see Sect. I4.1|l allows us to apply our method to compute 
benchmark solutions to test the accuracy of other methods to solve the coupled angular momentum and toroidal induction 
equations. A remarkable advantage of our approach is that it allows us to calculate the solution at arbitrarily spaced intervals 
of time, without the need to match Courant-like criteria, which is usually a severe limitation for other numerical methods of 
solving parabolic PDEs. 

The problem we are considering does not manifest the tendency to develop jump discontinuities, akin shocks in gas 
dynamics. This is another consequence of the parabolic nature of the equations we are solving, which tends to smooth out 
gradients with time. Therefore, our calculations are not significantly affected by the Gibbs phenomenon^ which can severely 
limit the quality of solutions computed by means of spectral methods. 



Applying a spectral method similar to ours to a problem already solved with different techniques, ICallvl (|199ll ) showed 
that the accuracy of the solution does not degrade with time as far as the number of eigenfunctions included is greater than 
the number of modes actually excited. With a resolution of TV = 50, K = 50, we are able to properly model the evolution 
of our system when we adopt an enhanced diffusivity factor R E s = 10 4 , corresponding to a Reynolds number of the order of 
10 8 for the solar core (note, however, that our code w orks under the as sumption of axisymmetry) . For comparison, the best 
numerical MHD codes introduced so far, such as ASH jciune et al.ll 19991) which is also based on a spectral method, can reach 
Reynolds numbers of the order of 10 5 — 10 6 in fully 3D calculations jZahn et al. 2007 ). 



4 RESULTS 

To illustrate the results of our modelling it is convenient to define a reference case, constructed with representative values of 
the various parameters. Specifically, we assume n = 0.1 Rq, = 0.7 Rq, with an initial rotation profile having fio = 25 Qq, 
Afi/fio = 0.01, r stC p = 0.5 Rq, Ar step = 0.2 Rq (cf. Eq. (|10|l ). A confined dipole configuration (see Eq. |[6|) is used for the 
poloidal field, with an amplitude Bo = 1 G. 
A result of previous works 

rotation of the radiative core within the solar age. Therefore, we assume an effective viscosity enhancement factor R c g = 10 4 , 
as defined in Eq. ([2|. It suffices to reproduce a rigid rotation in the core at the age of the Sun in our calculations, in agreement 
with helioseismic results. 

For the reference model, we adopt N — 50, K = 50, and a grid of 1000 points in the radial and meridional directions 
for the numerical evaluation of the radial and angular eigenfunctions and their scalar products, respectively. The number of 



CM93l ; lRK96h is that molecular viscosity alone is not capable of reproducing the nearly rigid 



1 The routines are available at http://www.netlib.org/lapack/. 

2 The partial sum of a Fourier series prone to the Gibbs phenomenon systematically overshoots the true value of the function with 
spurious oscillations in the neighbourhood of any jump di scontinuity. The problem is not alleviated by an increase of N and K because 
it is intrisic to the adopted representation of the solution llMorse fc F cshbach 1 9531) . 
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Figure 1. Upper panel: Evolution of the model without wind braking for resolutions N X K = 20 X 20 (dashed line), 30 X 30 (dotted 
line), and 50 X 50 (solid line), respectively. The extreme values of the energy E^t and £™ ln , as given by Eqs. J33D and 1341 . are also 
shown (dash-dotted lines). Lower panel: absolute difference between the reference resolution (50 X 50) and the 20 X 20 (dashed line) and 
30 X 30 (dotted line) solutions. 



grid points was varied in the range 500 — 1000 to check the consistency. The current choice of 1000 points was found to be 
sufficient for the quadrature routine to reach the accuracy of ~ 10~ J in relative units. The solution is calculated at 800 time 
ins tants, logarit hmically spaced from the ZAMS to the present solar age (for more details on the numerical implementation, 



Spadalboosh . 



Note that to obtain an accurate solution, it is necessary that the radial "eigenf unction resolution", i.e., (ra —ri)/K, be 
smaller than n. With n = 0.1 Rq, and the choice K ^ 30, this condition is safely realised. 



4.1 Conservation of energy and estimate of the truncation error 

Following Sect. 13.51 we define a global truncation error as: 

4k = Et - ^ ^ \ u nk + ( -r- ) Ak \ = 1 + ( 7^ ) r ■ 

n=0k = I ^ AJ J n = iVfc = K I ^ ' J 

From the last equality, we see that a% K is simply the sum of the squares of the coefficients of the neglected terms, and is thus 
always positive. 

Test calculations were performed without wind braking, an initial differential rotation AQ./Q.Q = 0.25, and N = K = 
20, 30 and 50, respectively. This suffices to test the code performance at different resolutions (i.e. with different values of N 
and K), focusing on the role of the coefficients corresponding to the smallest spatial scales. Testing the complete solution, 
namely including wind braking, is in principle feasible, but this would essentially influence only the woo term. 

As discussed in Sect. 13.51 in our high Reynolds number regime there is an initial phase, lasting ~ 3 • 10 5 yr, when total 
energy is conserved to a high degree (see Fig. [TJ. Then dissipation begins to decrease the energy down to the minimum value 
set by angular momentum conservation. A complete dissipation takes place within 10 s yr, independently of the adopted spatial 
resolution. 

Concerning the estimate of the truncation errors, all the calculations exhibit an excellent agreement (with unk ~ 10 -4 
in relative units) with the exact values of the initial and final energies and E™ ln during the initial and final stages of the 
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evolution, respectively. They can be computed from: 
and 



pr 



f2 in (r) 



n 



dr 



1 "-^core (0) 



_ 2 / 

with /core the total moment of inertia of the core scaled to 7q, i.e.: 



(33) 



(34) 



^corc — ~ 



pr 4 dr. 



The comparison of the less resolved runs with the reference run, 50 x 50, is particularly interesting during the intermediate 
phase, i.e. between ~ 10 6 — 10 7 yr, when the largest differences appear (see, in particular, the lower panel of Fig. [T]). Apparently, 
our reference solution reproduces more accurately the dissipation rate, given by the slope of the corresponding plot in the 
upper panel of Fig. [1] owing to the inclusion of a larger number of terms in the series. After about 2/3 of the total energy 
dissipation has taken place, a change of the slope occurs also in the reference solution. If this is a numerical artefact due 
to the truncation of the energy transfer to smaller spatial scales, then our computed time-scales for relaxation to a rigid 
rotation state could be slightly overestimated. However, this is not a serious limitation if we are in terested in an order of 
magnitude est imate of the coupling time-scale akin to, for example, the parameter r c in the model by Allainl 1 1998h (see also 
Bouvier 20081 ) . A suitable correction factor could, in principle, be estimated from Fig. [TJ by extrapolating the slope of the 
energy variation vs. time beyond ~ 2 • 10 6 yr. With this simple approach, we can estimate that the e-folding time for energy 
dissipation, as derived from the upper part of the plot of the energy vs. time, is overestimated by no more than a factor of 
2-3. 



4.2 Evolution of the reference model 



The model with wind braking evolves through the following three main phases: 
a) Linear build up of toroidal field 

The toroidal field amplitude, starting from zero with the initial condition in Eq. grows with time for a few tn, with 



tn <C tA ~ 3.5 • l(r yr for B 



1 G. In this phase, advection dominates over diffusion in the induction equation (|15[) . so that: 
dB v 



dt 



r sini9B„ ■ Vfi. 



Noting that, in the initial phases, Vfi = 4^r (cf. Eq. (JTOjl) and integrating with respect to the time: 



QQ 

r sin fiB r — At. 
or 



(35) 



The stretching of the poloidal field lines by differential rotation is characterised by the smallest time-scale, namely tn- The 
amplitude of B v increases linearly with time, feeding on the shearing of B p , and is proportional to the initial differential 
rotation Att and poloidal field amplitude Bo- The toroidal field is antisymmetric with respect to the equator, as shown in the 
contour plot of Fig. [2] This is a consequence of Eq. (|35p . since Q is initially independent of i? and B r ~ siniJcosif, as can be 
seen from Eq. ((6} for "J. 

The linear growth of the field eventually ends when a balance is established between the torques due to wind braking and 
the azimuthal component of the Lorenz force, which reacts to further shearing, 
b) Torsional Alfven waves 

The second phase, ranging approximately from 10 3 to 10 6 yr for Bo = 1 G, is characterised by the excitation and progressive 
damping of torsional Alfven waves. When the Lorentz force becomes effective, it acts as a restoring force on any further 
stretching of the poloidal field lines. This excites oscillations of B v and fl, propagating along each poloidal field line (cf. Sect. 
13.41 and Fig. [3}. These waves have a linear behaviour for arbitrary amplitude, because their phase velocity, Va p = Bp/^/inp, 
is independent of the toroidal field intensity or the rotation rate. Neighbouring poloidal field lines, however, oscillate with 
slightly different phase velocities because |B P | and p are not uniform in the core. This leads to the development and progressive 
increase of a phase lag between waves propagating along neigh bour field lines. Eventual l y these wave s get in opposition of 
phase and dissipate quickly. This process is called phase mixing jSpruitlll987l : ICallvlll99ll : [Spruitlll999h . 

Some interesting features of the oscillations are illustrated in Fig. [4] Control points located at different depths show that 
oscillations in the field amplitude are more readily excited and damped in the external part of the core, i.e. closer to the 
driving torque (Fig. 4(a) I. The differential rotation imposed as an initial condition is significantly smoothed out before the 



wind braking becomes important (Fig. |4(c")| . Control points located at the same radius were chosen to investigate the effect 
of a trapping of the oscillations inside the so-called "dead zone" of the poloidal field, i.e. in the domain near the magnetic 
neutral "O" point at r ~ 0.4 Rq on the equator. At lower latitudes, the B v oscillations persist for a longer time (Fig. 4(b) l 
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Figure 2. Isocontour plots of O (left panel, in units of Qq) and B v (right panel, in G) after 10 3 t\. 
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Figure 3. As in Fig. [2] for t = 3.1 ■ 10 3 yr ~ 1 t\. Out-of-phasc oscillations on neighbour poloidal field isosurfaces are clearly evident 
in the B v isocontours in the right panel. 
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(b) B v at fixed radius. 
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Figure 4. B v (upper panels) and Q (lower panels) vs. time at some control points. On the left panels, points at r = 0.2, 0.4, 0.6 -Rq, 
with 1? = 7r/4 fixed, are indicated by dotted, solid and dashed lines, respectively. On the right panels, points at $ = n/12, n/A, 57r/12, 
with r = 0.4 Rq fixed, are indicated by dotted, solid and dashed lines, respectively. Dotted and dashed curves in |4(a)| and |4(b)| were 
given an offset of ±100 kG for the sake of clarity. 



and rotational braking is delayed (Fig. 4(d) I. We conclude that, with this choice of R c s, oscillations are already damped in a 
considerable fraction of the core, including the dead zone, before the wind braking of the star as a whole becomes important. 

Our analytic treatment allows us to evaluate precise upper and lower bounds for the total energy contained in the subgrid 
lenghtscales during the phase mixing process. For instance, considering the cases plotted in Fig. \T\ and that the total energy 
is always greater than the sum of the squared amplitudes of the included modes, we find that 2.4 x 10 -3 ^ aff K (t)/ET(t) < 
9.1 x 10" 3 at any time t during the phase mixing evolution. This rigorous treatment of the subgrid contribution is a new and 
interesting feature of our analytic approach, 
c) Quasi- stationary evolution 

For the rest of the computation, the evolution proceeds at a slower pace, in a quasi-stationary regime reminiscent of Ferraro 
isorotation, i.e., with angular velocity almost uniform along each poloidal field line. Actually, this is the expected outcome of 
the phase-mixing of torsional waves, which smooths out the fluctuations of angular velocity on poloidal field isosurfaces, that 
is the surfaces where ^ = const. Neglecting the diffusion processes, the stationary condition is reached when: 



Bp • Vfi =0 
Bp • V(rsint?B v ) = 0. 



(36) 



Eqs. (|36[) predict th at the rotation rate must b e constant on ^ isosurfaces, and that the azimuthal component of the Lorentz 
force must be zero jGaraud fc Guervillvlliooi ). These conditions are achieved at late times in our model, as can be seen in 
Fig. [5] This is a consequence of the ratio of the time-scales involved, because deviations from the conditions in Eqs. (136[1 are 



produced on the wind braking or diffusi ve time -scale, i.e., from ~ 10 to ~ 10 yr, but are compensated for on the much 



shorter AlfVen time-scale £a ~ 10 3 yr (cf. CM93])- Nevertheless, the magnetic (and viscous) stresses are still acting to transfer 



angular momentum on the long time-scale characteristic of the late wind braking process. 

As already noted, a model with R c g = 1 does not attain a condition of uniform rotation in the core within the solar age, 
even in the presence of a large scale poloidal field. Magnetic transport of angular momentum due to phase mixing is effective 
only along poloidal field lines, while the only means to couple the plasma across magnetic surfaces is by the action of the 
viscosity. Molecular viscosity alone fails to ensure an effective rotational coupling within the dead zone of the poloidal field 
because its coupling time-scale is of the order of R 2 /v, i.e., longer than the age of the Sun. In the model considered here, it 
is an enhanced viscosity which eventually enforces rigid rotation, as seen in Fig. [6] 
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Figure 5. As in Fig. [2] for t = 1.23 Gyr or 350 000 tA- H (left panel) and By, (right panel) isolines do not change for the rest of the 
evolution, apart from an overall scale factor. Note that £7 has become nearly constant on poloidal field lines and B v = inside the "dead 
zone", where rsm-DB^ must be both antisymmetric with respect to the equator and constant along field lines. 
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Figure 6. Upper panel: latitudinal average of the angular velocity vs. time at three different depths: T2 (solid line), r\ (dashed line), 
and r c = 0.4 Rq (dotted line), in the vicinity of the neutral magnetic point. Lower panel: relative difference between the average angular 
velocities at and r c vs. time. 
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Figure 7. Angular velocity as a function of the fractional radius r, at the equator (solid line) and at 45° latitude (dashed line), at the 
age of the Sun (~ 4.6 Gyr), for R c ff = 10 (left panel) and R e g = 10 4 (right panel). 



Note that the ev olution of the average angular velocity is very similar to that obtained by RK96I (see the upper panel of 
Fig. [6] and Fig. 4 of IRK96I) . apart from a shallow dip in our curve around 10 6 yr due to the inclusion of a saturation for the 
angular momentum loss in the wind braking law adopted in our model. 

The two panels in Fig. [7] show plots of the angular velocity as a function of the fractional radius for two representative 
latitudes, corresponding to viscosity enhancements R c g = 10 and R c g — 10 4 , respectively. Fo r comp arison purposes, they are 
presented with the same axes of the upper panel in Fig. 6 and the central panel in Fig. 7 of RK96I . respectively. The profile 
near the pole is not shown here because the exclusion of the central part of the core (r < n) produces spurious oscillations 
near the axis that are not damped effectively with a viscosity value as low as R a s. Given the negligible amount of angular 
momentum near the axis, however, this does not influence the whole calculation significantly. 

The error on the angular velocity was estimated for both cases according to the method of Sect. |4~T1 obtaining a relative 
error < 4% for R ctt = 10 (corresponding to a Reynolds number Re ~ 10 11 ) and < 2.5% for R aH = 10 (Re ~ 10 s ). The 
appropriate error bars are shown in Fig. [7] 

The small oscillations in the left panel of Fig. [7] are well within the error bars and are a consequence of the difficulty to 
represent a nearly co nstant plateau with a finite number of periodic eigenfunctions. Overall, both panels compare well with 
their counterparts in IRK96I . apart from the immediate vicinities of the radial boundaries n and ri, where the effect of the 
different boundary conditions is more marked. The stronger coupling between the equator and 45° that can be seen in our 
results is a consequence of the different prescription for the poloidal field flux function. 

If conditions close to Eqs. (|36[1 are eventually established, a poloidal field having a multipole configuration of an order 
higher than dipolar could be mor e efficient in imposing rigid rotation in the whole stellar core because its dead zones are 
significantly smaller 1 Spruit 1999). This is in fact observed in the case of our quadrupolar model, that attains a solid body 
rotation within ~ 1 Gyr (see Se ct. 14. 3[1 . Such a time-scale is to be compared with a decay time of the quadrupolar mode of 
2.7 Gyr, as estimated by lRK9fj . 



4.3 Differential rotation of the core 

To ease the comparison with previous works, a measure of core differential rotation is defined as: 

D(t) . f * f * W'W-^'Wdrdca**. (37) 

V ; 2{rl~r{) J_, J ri fi(r 2 ,7r/2,t) V ' 

In Fig. [8] we show the evolution of D(t) for some models compared with the reference model. 

The curves corresponding to R e g — 10 2 — 10 3 show that such low values of the viscosity are not enough to attain a 
uniform rotation within the age of the Sun. Notably, decreasing R c g has the effect of attaining the maximum level of internal 
decoupling at later times. 

The poloidal field geometry has a remarkable impact on the phase-mixing phase and on the quasi-stationary regime 
that is subsequently established. Apart from the changes in the shape of the angular velocity and toroidal field isolines, 
an assigned flux function with quadrupolar symmetry determines both qualitative and quantitative modifications in the 
differential rotation evolution with respect to that with a dipolar symmetry. As it is shown in Fig. [8] the decoupling phase 
is shortened and D(t) presents a plateau around 10 6 — 10 7 yr in the quadrupolar model. Conversely, the role played by the 
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Figure 8. Differential rotation measure D vs. time for different models. Solid line: reference model (dipolar poloidal field with Bo = 1 G 
and R B g = 10 4 ); dashed line: model with R c g = 10 2 ; dash-dotted line: model with i? c ff = 10 3 ; long-dashed line: model with i? c ff = 10 s ; 
dash-triplc-dotted line: model with quadrupolar poloidal field (Bo = 1 G) and R e ff = 10 4 ; dotted line: non-magnetic model (i.e., Bq = 0) 
with i? c ff = 10 4 . The oscillations of D near 10 4 and 10 7 yr are due to torsional Alfven waves. 



poloidal field intensity Bo is quite modest, as is apparent in Fig. |9j because the winding up by differential rotation leads in 
any case to similar azimuthal stresses B r B v at the end of the linear build u p phase of the azimuthal field. 

Our reference model behaves very similarly to that computed by CM93I with similar parameters for ages > 5 • 10 7 yr. The 
difference at earlier ages is due to the fact that we assume a differential rotation on the ZAMS. 

Varying Afi/flo in Eq. (|10[) in the range 0.01 to 0.1, we verified that the amplitude of differential rotation after ~ 10 7 yr 
does not change significantly. This is consistent with a weak dependence on the details of the rotation profile on the ZAMS 
at later stages. 

Our results could be particularly interesting to explain the differences between fast and s low rotators in the core-envelope 
coupling time-scale, according to the phenomenological models by, e.g., Allain ( 19981 ) and Bouviei ( 20081 ) . If fast rotators 
reach the MS with a predominatly quadrupolar field, then their coupling time is about ten times shorter than that of slow 
rotators, that we may assume to have a predominantly dipolar fie ld. The selection of the preferred mode may be determined 
by the stellar dynamo during the PMS phase (see, e.g., Sect. 3 of Moss et al. 20081 ). 



5 DISCUSSION 



The method used here shares its mathematical foundations with that applied by lLanzal 1)20071 ). where the solution of the 
PDE expressing the conservation of the angular momentum was expanded in a generalised Fourier series. In this work, we 
further develop this approach to solve two coupled PDEs, with non-homogeneous boundary conditions. This method offers 
several advantages over finite difference methods, notably the fact that the accuracy does not degrade with time and that the 
calculation of the solution at ev ery previous step is not needed if only one particular instant is required. 

According to Spruit ( 19991 ). the strongest instability in a magnet ized stellar radiat ive zone should be that described 
bv lPitts fc Tavleil (|l985l ). Spruit's analysis was numerically tested by IZahn et al.l ((2.007), using a 3D MHD spectral code. 
They confirmed that the strongest non-axisymmetric unstable mode has an azimuthal wavenumber m — 1 as for Pitts- Tayler 
instability, but there was no hint of a regenera tion of the poloidal magnetic field owing to a dynamo action associated with 
the instability, as conjectured by ISpruiti (120021 ). Non-axisymmetric motions associated with the instability seem to behave 
as Alfven waves rather than turbulence. Their contribution to the transport o f ang ular momentum appears to be negligible, 
contrary to the conjectures of Spruit ( 19991 ) , Denissenkov fc Pinsonneault ( 2007 ) , and Denissenkov et al. ( 20081 ) , who envisaged 
an increase of the effective turbulent viscosity and diffusivity as a consequence of the motions associated with the instability. 
Therefore, our assumption of a purely axisymmetric model to study angular momentum transport within a radiative core 
gains support. 

The choice of the initial conditions used in our calculation was mainly motivated by simplicity reasons because no 
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Figure 9. D vs. time for different poloidal field intensity. Solid line: reference model (Bo = 1 G); dotted: So = 0.1 G; dashed: Bo = 0.01 
G; dash-dotted: non-magnetic model. The oscillations of D near 10 4 and 10 7 yr are due to torsional Alfven waves. 



observational constraints on the radial rotational profile or inner magnetic field configuration of stars on the ZAMS is presently 
available. Nevertheless, except for the geometry of the poloidal field, our choice of the initial conditions does not significantly 
affect the solution and it is, therefore, not critical for the angular momentum transport problem we have addressed. 

Neglecting the poloidal magnetic field diffusion is acceptable on the basis of the quite long time-scale of this process. This 
could be questionable if some source of enhanced m agnetic diffusivity is present in the core, perhaps in association with the 
still elusive viscosity enhancement, as sugg ested by|lK96. The effect of the diffusion of the poloidal field is merely a decay 
of its amplitude. This has little influence on the solution and was not taken into account to avoid unnecessary complications. 
On the other hand, a possible spatial reshaping of the poloidal field lines, could have a greater impact on the phase-mixing 
process and the late stages. However, the inclusion of such a field evolution would require a knowledge of the initial magnetic 
field configuration, which is not available. 

Finally we note that a poloidal field in the core may imprint the differential rotation at the outer boundary of the core 
into the radiative interi or of the Sun, contr ibuting to a significant downward spreading of the tachocline which is contrary to 
helioseismic results (cf. Brun fc Z ahn 2006). Neglecting the diffusion of the poloidal field and assuming a latitude-independent 
radial gradient of the angular velocity at the outer boundary, as in our calculations, prevent such an effect. This is equivalent to 
assume that the tachocline is not dynamically connected with the core poloidal field, which could be justified if the tachocline 
is confined into a thin layer above r — T2 by other effects, e.g., a highly anisotropic turbulent diffusivity, as proposed by 
Spiegel fc Zaiml l|l992l ). 



6 CONCLUSIONS 

We introduce an exact analytic solution of the coupled equations of angular momentum and azimuthal magnetic field evolution 
in a radiative core under the main assumptions of axisymmetry and negligible meridional circulation. From such a solution 
we derive a numerical spectral method that allows us to address the problem of angular momentum transport in stellar 
radiative regions. Our illustrative solution of the angular momentum evolution in the radiative core of the Sun is compared 
to previous numerical models to discuss the advantages of the present approach. It is particularly intesting because it allows 
us to define an analytic upper bound for the numerical errors and it can easily reach hydromagnetic regimes characterized 
by Reynolds numbers of the order of 10 7 — 10 s , that are not accessible with other numerical techniques. Our approach 
is particularly interesting because it provides a rigorous treatment of the kinetic and magnetic energy distribution among 
different lenghtscales during the phase mixing process. Previous numerical methods were not capable of a rigorous evaluation 
of the contribution of the subgrid spatial scales, while our method does it by applying the mathematical closure expressions 
for generalized Fourier series and the conservation of energy. 

With our approach, we studied the evolution of the rotational decoupling of the core, defined through a measure of its 
differential rotation. Our results confirm that the uniform rotation of the core of the present Sun, as deduced from helioseismic 
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inversions, can be reproduced only if an enhanced viscosity, ~ 10 4 times greater than the molecular value, is assumed. This is 
an early finding of previous works, but the identification of the underlying physical process(es) is still debated. Asteroseismic 
measures are needed to infer how general the rigid rotation of the core of the Sun is. 

Another interesting result is that a quadrupolar poloidal field leads to a coupling time-scale about one order of magnitude 
shorter than a dipolar field. We conjecture that such a difference in the geometry of the poloidal field may be a consequence of 
the different dynamo regimes operating in fast and slowly rotating stars during the PMS stage, respectively. If verified by future 
studies, this may explain the different coupling time-scales required by the phenomenological models of rotation braking of fast 
and slow rotators, respectively. We intend to address further such an interesting topic in a forthcoming work. Moreover, our 
analytic treatment of the coupled equations of the angular momentum and magnetic field evolution in an axisymmetric MHD 
system can be applied to other astrophysically relevant problems, e.g., to study the torsional oscilla tions of a magnetized shell . 
Thi s problem has interesting applicatio ns to, e.g., the oscillations in roAp stars, as discussed by, e.g., Rincon fe Rieutord ( 20031 ) 
and Rees e. Rincon. fe Rieutord! (|200J). They have considered the case of a magnetized shell of uniform and constant density 
bathed by a dipolar magnetic field in the absence of rotation, studying both the poloidal and the toroidal axisymmetric modes 
as well as non-axisymmetric modes by means of linearized equations, i.e., considering only small perturbations. Our approach 
allows us to study only torsional Alfven waves, i.e., the toroidal modes in their classification, but without the simplifying 
assumption of small amplitude, and considering also the stratification of the medium and stellar rotation, when the axis of 
the dipole field is aligned with the rotation axis. Appl ications of our approach to strongly magnetized, degenerate stars, suc h 
as magnetic white dwarfs or neutron stars (cf., e.g., lOkita fe Koiimall2005l : iGlampedakis. Samuelsson. fe Anderssonlbood ). 
may also be of interest. Other possible applications include an implementation for cylindrically symmetric systems, e.g., an 
accretion disc around a protostar threaded by a large-scale dipolar field, to study the torsional Alfven waves and their phase 
mixing within the ionized region of the disc. However, in this case our model represents an even stronger idealization than 
in the cases mentioned above because the geometry of the poloid al field and the large-scale poloidal f low shown by detailed 



numerical simulations are far from being simply dipolar (cf., e.g., von Rekowski fe Brandenburd 2004 ) 
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APPENDIX A: THE WIND BRAKING BOUNDARY CONDITION 



The amount of shear at the lower boundary of the CZ, assigned through the function W(t) in Eq. ([5J, can be determined 
using the integral form of the angular momentum equation applied to the volume of the core c €: 



d_ 

dt 



/pr sin ' mdr = / V ■ (pr sin iVVfi)dr 
J is 



2 • 2 

pvr sin 



dr 



dy, 



where ^ is the boundary of the core, i.e., the sphere of radius r 2 , and we have applied Gauss theorem, neglecting the azimuthal 
magnetic stress term B r B v because B v = at r = r 2 (cf. Eq. Q). 

The l.h.s. is the total angular momentum loss suffered by the core. Assuming that ^ on 5? does not depend on •& and 
performing the surface integration, we obtain: 



-p(V 2 )^(r 2 )r 2 



dr 



dJ 



Kawalei 1 1988h proposed an expression for the angular momentum loss by a magnetised wind, slightly modified by Chabover et al. 
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1 1995al lbh ~,o account for the observed saturation at high rotation rates: 



dt) w - Kw \m q ) \r q ) n,>n c > (AL) 

where K w ~ 10 47 g cm 2 s is a constant calibrated in such a way to obtain the present surface angular velocity of the Sun. 
The adopted braking law includes a basic dependence on the surface rotation rate fi* by simple power laws, along with a 
dependence on stellar parameters 7?*, M*. 

To assign the boundary condition in the second of Eqs. JS), we equate the flux of angular momentum outside of the core 
to the angular momentum lost in the wind as given by Eq. (|Alj) : 

dt) c ~\dt 

In conclusion, we have in the solar case (M, = Mq, R, = Rq) 

dQ, 
dr 



where /(f2*) is the law in the rightmost factor in the r.h.s. of Eq. (| Al fl . 

The dependence on fi*(t) introduces a non-linearity in the problem. We use the equatorial value of the angular velocity 
fi(r2, § ,~t) to represent f2,(t) for the computation of the boundary condition. Note that we are implicitly assuming that all 
the angular momentum lost in the wind is extracted from the boundary of the core. This is equivalent to assume that the 
CZ and the tachocline can extract angular momentum from the core on a time-scale much shorter than the timescale for the 
angular momentum redistribution within the core itself. Given the high value of the turbulent viscosity in the CZ this appears 
to be a plausible hypothesis. 



APPENDIX B: SELF-ADJOINT OPERATORS AND EIGENFUNCTION EXPANSION 

Let us consider a linear differential operator of the form: 

C^] = V ■ (pV0) + q 0, (Bl) 

with p, q given functions of r and ■&. 

The definition of self-adjointness commonly applied to matrices can be generalised for this kind of operators. We introduce 
the Green formula for the operator C, i.e.: 

j {uC[v\ - vC[u\) d 3 r = j> p (uVv - wVw) ■ ndS, (B2) 

where u and v are arbitrary differentiable functions and the surface integration extends over the boundary S of the integration 
domain C, with n being the unit vector in the direction of the ourward normal to 5. It is easily seen that the operator C, as 
defined by Eq. l|Bl]l. satisfies Eq. (|B2|) . 

Let us consider two arbitrary differentiable functions (/>, tjj satisfying linear boundary conditions on S of the kind: 

c^ + /3V0-n =0, 

aip + pX/ip ■ n =0, 1 ' 

where a and /3 are arbitrary real numbers. Then, making use of Eq. ()B2|I with u — (j> and v = ifr, we find that: 

which is the defining property of a self-adjoint differential operator. 

The eigenvalue problem for a self-adjoint linear operator is specified by the equation: 

C\<j>n] + KnW(f>n = 

together with the boundary conditions in Eq. (|B3|) . It is also called a Stiirm-Liouville boundary value problem. For this kind 
of problem, an infinite set of eigenfunctions can be proved to exist, forming a complete and orthogonal set. 

The self-adjoint operators appearing in Eqs. (|16[) and (|17[1 can be put in the form of Eq. (|B1|) with the positions: 

p = pr 2 sin 2 

q = 0, 

for the Z problem, and 
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for the S problem. 

We shall now apply the Green formula to derive Eqs. (|26|l . With u = Q and v = Z, Eq. (|B2|) becomes: 



nV ■ (pr sin -3 MZ^) - Z nk V ■ (pr sin 13 *V!1) (f r = * pr v sin ■& fi 



> dZ nk 
dr 



Znk-^r- ) dS^ , 
or 



f'e J s? 

where the integration is extended over the computational domain V = [ri,T2] x [— X, +1] X [0, 2n], bounded by the surface 
S' = Si U S2, where Si and 5*2 are the two spherical surfaces of centre O and radius n, r2, respectively. This formula allows 
us to take into account the boundary conditions. Using Eqs. (I22p and the r.h.s. becomes: 

-1 



pr i^sin ■& I 



•5" 



dr 



" dr 



dT, 



-2n p(r 2 )u(r 2 )rU n k(r 2 )W(t) 



sin^ ■&P^ 1 ' 1 'd cos 1? 



-^nfc(i). 



The left hand side of the Green formula can be further manipulated by using Eqs. (|16[l and (I14|l . i.e. with the substitutions: 



V ■ (pr sin iWZ n fc) 
V ■ (pr 2 sin 2 tfi/Vfi) 



-pr sin tf\ nk Z nk , 

pr 2 sin 2 $-7^- — Bp • V(r sin #-B^ 



Recalling that Z nk (r,$) = ( nh (r)P n ' >($), we obtain: 



f OV • (pr 2 sin 2 iVVZ nfc )d 3 r = 
/ Z„ fc V ■ (pr 2 sin 2 iVVfi)d 3 r = 



pr 2 sin 2 '&C, nk P n 1 ' 1 >Q,r 2 dr dcosi? = ~\ nk u nk 



1 

1 r r 2 



-1 J r 
r + 1 r r 2 



2 ■ 2 QA nfl.l) 9f2 2 , , q 

pr sm vCnk"n -^-f dr acosw 
at 



/ pr 2 sin 2 tffnfcPn'^Bp • V(r sin , &B v )r 2 dr dcos-d 

/ Snkmhftmh J ■ 



dt 

The last two equalities were established using the definition of the generalised Fourier coefficients in Eqs. (|24p and (|25p and 
their time derivatives. Moreover, the linear dependence of the source term on B v allows us to express the spatial dependence 
of its expansion through suitable coefficients S„ km h, entering as factors of /3 m h(t) and whose explicit expression is given 
below. Thus, by applying the Green formula to Eq. (|16[) . provided that the eigenfunctions Z nk are conveniently normalised 
(i.e. {Znk, Z nk ) = 1), we obtain the first of Eqs. 



duj n l 

dt 



/ Snkmhftmh J — ^nfc • 
mh / 



Applying the same line of reasoning, the derivation of the second of Eqs. (|26[) is quite straightforward. The boundary terms 
vanish completely in view of Eqs. (|9} and p2[) and we are left with: 



with the shorthand notation: T>[z. n k] = V • {r]X7~„k 



{B v V[~. n k\ - E n kV[B v ]} <f r = 0, 
rf 1 



r 2 sin 2 § 

Substitution of the eigenvalue problem, Eq. (|17[1 . and of the PDE, Eq. (|15[) , i.e., 



T>[~ nk ] 
V\B V ] 

leads to: 



dB^ 



dt 



r sin #Vfi 



r r 2 1 2 

Pnfc J J £,nkP n B v r dr dcos-d 

/ r + l rr 2 Qr> r + 1 rr 2 

'^ V \J J ^Pn-^r 2 dr dcos$-J J £ n kP„r sin $B P ■ X7ttr 2 dr dcos ■& 



-/J-nkPnk - M 



, / df5 nk _ ^ 

" I dt ^ n 

\ mh 



kmh^mh — 0, 



which is indeed equivalent to the second of Eqs. (1261) . 

For reference, we provide here the complete expressions of S„kmh and T n kmh, omitting the huge amount of algebra 
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necessary to derive them 
hl d 



y/l-xt-g-pWpldx ■ I <£— [rUh] tnkdr 



' nkmh 



dx 
-dP (1 ' 1} 



rip— 3 — t,nkdr 
dr 



where the functions (p and x appear in the expression of the poloidal field flux function, i.e., \l/ = tp(r)x{x), with x = cost?. 

From these expressions, a remarkable property of the matrix A in Eq. (|27[) can be proved. Using the boundary conditions 
for the eigenfunctions £ n fc in Eqs. (|22[) and integrating by parts, we find: 

Tnkmh = —Smhnk- (B4) 

With the notation L = diag(— \j/M v ), M = diag(— fij/.'M^), § = (S^;), where j, / are indexes that number the orderer couples 
(n, fc), (m, h), respectively, e.g., j — Kn + k, I — Km + h, we can write: 

+ 





s 




M 



L 







M 






§ 








Thus we see that A is the sum of a diagonal matrix plus an antisymmetric matrix. Since both diagonal and antisymmetric 
matrices commute with each other and with their respective transposes, i.e., they are normal matrices, also their sum A is 
a normal matrix. This is a very useful property, because for normal matrices a generalisation of the spectral theorem holds, 
which is then applicable to Eq. (|27[) . as discussed in Sect. 13.41 



